R Notebook

library(tidyverse)
library(mgnet)

mg <- readRDS("cache/gs_netw_filt.rds") %>%
  filter_taxa(!(comm_id %in% "7"))
mg

==== mgnet Object Summary ====
General Info:
  Samples: 912
  Taxa: 1235
  Only normalized data available.
  Sample Meta Info: Region, country_alt, city, year , etc... 
  Taxa Meta Info: source, type, kingdom_ncbi, phylum_ncbi , etc... 
  Network: 26071 edges, ~3.42% density
  Detected Communities: 6
  Community Sizes: 318, 247, 126, 113, 101 
  Isolated Nodes: 0 
==== End of mgnet Object Summary ====

# Create visual properties of the network, one by one
mg <- mg %>%
  # Color at family level for mOTU
  mutate_taxa(family_fct = if_else(is.na(family_ncbi), 
                                   "Unknown", family_ncbi)) %>%
  mutate_taxa(family_fct = fct_lump_n(family_fct, 20)) %>%
  colorize_taxa(family_fct, color_to = "color_family", 
                colorspace = list(h = c(0, 360), s = c(.2, .6), l = c(.7, .9))) %>% 
  mutate_taxa(color_family = case_when(
    family_fct == "Unknown" ~ "#D3D3D3",
    family_fct == "Other" ~ "#F7F7F7",
    TRUE ~ color_family
  )) %>%
  # I switch manually at two important families to make it more evident in the 
  # final results
  mutate_taxa(color_family = case_when(
    family_fct == "Pseudomonadaceae" ~ "#8FA4DAFF",  
    family_fct == "Neisseriaceae" ~ "#EBE7DAFF",
    family_fct == "Enterobacteriaceae" ~ "#E393D9FF",
    family_fct == "Aeromonadaceae" ~ "#CBA8B7FF",
    TRUE ~ color_family
  )) %>%
  # Color for amr classes
  mutate_taxa(class_amr_fct = factor(class_amr)) %>% 
  mutate_taxa(class_amr_fct = fct_lump_min(class_amr_fct, 10, other_level = "['other']"), .by = "") %>%
  colorize_taxa(class_amr_fct, color_to = "color_class_amr") %>%
  mutate_taxa(color_class_amr = if_else(class_amr_fct == "['other']", "#F7F7F7", color_class_amr)) %>%
  # Create a new column color_nodes where summarize amr and mOTU colors
  mutate_taxa(color_nodes = if_else(source == "mOTU", color_family, color_class_amr),
              color_nodes = if_else(source == "mOTU", 
                                    adjustcolor(color_family, alpha.f = .85),
                                    color_class_amr)) %>%
  # Color for the edges of the nodes
  mutate_taxa(vertex.frame.color = case_when(
    type == "mOTU" ~ "#838B83",
    type == "acquired" ~ "#8B0000",
    type == "functional" ~ "#008B45"
  )) %>%
  # Communities colors
  mutate_taxa(comm_id_fct = factor(comm_id)) %>%
  colorize_taxa(comm_id_fct, color_to = "color_comm", 
                colorspace = list(h = c(0, 360), s = c(.5, 1), l = c(.4,.7))) %>%
  mutate_taxa(color_comm = adjustcolor(color_comm, alpha.f = .6))

# Create named vector for the colors useful to create the legends
color_family <- gather_taxa(mg) %>%
  select(family_fct, color_family) %>%
  distinct() %>%
  deframe() 

color_amr <- gather_taxa(mg) %>%
  select(class_amr_fct, color_class_amr) %>%
  filter(!is.na(class_amr_fct)) %>%
  distinct() %>%
  deframe()
names(color_amr) <- names(color_amr) %>%
  str_remove_all("\\[") %>% str_remove_all("\\]") %>% str_remove_all("'")
#color_amr["other"] <- "#F7F7F7"

color_comm <- gather_taxa(mg) %>%
  select(comm_id_fct, color_comm) %>%
  filter(!is.na(comm_id_fct)) %>%
  distinct() %>%
  deframe()
color_comm <- color_comm[as.character(1:ncomm(mg))]

# generic star vertex shape, with a parameter for number of rays
# respect the igraph example i add the frame color and width
mystar <- function(coords, v=NULL, params) {
  vertex.color <- params("vertex", "color")
  if (length(vertex.color) != 1 && !is.null(v)) {
    vertex.color <- vertex.color[v]
  }
  vertex.size  <- 1/200 * params("vertex", "size")
  if (length(vertex.size) != 1 && !is.null(v)) {
    vertex.size <- vertex.size[v]
  }
  vertex.frame.color <- params("vertex", "frame.color")
  if (length(vertex.frame.color) != 1 && !is.null(v)) {
    vertex.frame.color <- vertex.frame.color[v]
  }
  vertex.frame.width <- params("vertex", "frame.width")
  if (length(vertex.frame.width) != 1 && !is.null(v)) {
    vertex.frame.width <- vertex.frame.width[v]
  }
  norays <- params("vertex", "norays")
  if (length(norays) != 1 && !is.null(v)) {
    norays <- norays[v]
  }

  mapply(coords[,1], coords[,2], vertex.color, 
         vertex.frame.color, vertex.frame.width,
         vertex.size, norays,
         FUN=function(x, y, bg, fg, lwd, size, nor) {
           symbols(x=x, y=y, bg=bg, fg=fg, lwd=lwd,
                   stars=matrix(c(size,size/2), nrow=1, ncol=nor*2),
                   add=TRUE, inches=FALSE)
         })
}
# no clipping, edges will be below the vertices anyway
add_shape("star", clip=shape_noclip,
                 plot=mystar, parameters=list(vertex.norays=5))

Start Plots one by one

Community plot network

x <- mg %>%
  select_link(weight > 0) %>%
  group_taxa(comm_id) %>%
  mutate_netw(eig_pos_comm = eigen_centrality(netw)$vector) %>%
  filter_taxa(eig_pos_comm > quantile(eig_pos_comm, 0.25)) %>%
  mutate_netw(deg_comm = degree(netw)) %>%
  filter_taxa(deg_comm >= 3) %>%
  deselect_link() %>%
  ungroup_taxa()


set.seed(42)
layout_x <- layout_signed(netw(x), layout = with_graphopt(
  spring.length = 1,
  mass = 100,
  charge = .003
))

x_inv <- x[, ntaxa(x):1]
layout_x_inv <- layout_x[ntaxa(x):1, ]

png(filename = "intermediate/network_community.png", 
    width = 4200, height = 4200, res = 300)

par(mar=c(0,0,0,0))
plot(x,
     layout = layout_x,
     vertex.label = NA,
     posCol = rgb(0, 0, 1, .2), negCol = rgb(1, 0, 0, .2),
     widthFactor = .75,
     vertex.shape = if_else(pull_taxa(x,source) == "mOTU", "circle", "star"),
     vertex.size = colMeans(norm(x)),
     maxSize = 6, expFactor = 0.5, sumConst = 1,
     vertex.frame.width = 3,
     vertex.color = if_else(pull_taxa(x,is_from_gut), "#CD6839", adjustcolor("#FFFACD", alpha.f = .5)),
     vertex.frame.color = pull_taxa(x, color_comm))

legend(
  x = 1.1, y = -.92, xjust = 1,
  legend = names(color_comm),
  col = color_comm,
  pch = 1, cex = 1.9,
  pt.cex = 2.5, bty = "n", lwd = 3,
  title = "Communities", horiz = T
)

legend(
  x = -1.1, y = -.92, xjust = 0,
  legend = c("Human mOTU", "Other mOTU", "amr"),
  col = c("#CD6839",  adjustcolor("#FFFACD", alpha.f = .5), "black"),
  pch = c(19,19,11), cex = 1.9,
  pt.cex = 2.5, bty = "n",
  title = "Source", horiz = T
)
dev.off()

null device 
          1 

png(filename = "intermediate/network_taxonomy_only.png", width = 4200, height = 4200, res = 300)

par(mar=c(0,0,0,0))
plot(x,
     layout = layout_x,
     vertex.label = NA,
     posCol = rgb(0, 0, 1, .2), negCol = rgb(1, 0, 0, .2),
     widthFactor = .75,
     vertex.shape = if_else(pull_taxa(x,source) == "mOTU", "circle", "star"),
     vertex.size = colMeans(norm(x)),
     maxSize = 6, expFactor = 0.5, sumConst = 1,
     vertex.color = pull_taxa(x, color_nodes),
     vertex.frame.color = pull_taxa(x, vertex.frame.color),
     vertex.frame.width = if_else(pull_taxa(x,source) == "mOTU", .5, 2))



dev.off()

null device 
          1 

Plot with only the legends of the networks

png(filename = "intermediate/netw_legend.png", width = 5250, height = 2100, res = 300)
plot.new()
par(mar = c(0,0,0,0))

legend(x = -.082, y = 1.1, xjust = 0,
       legend = tolower(names(color_family)), fill = color_family, title = "mOTU family", 
       y.intersp = 1.2,
       ncol = 2, bty = "n", cex = 2, title.adj = 0)

legend(x = .55, y = 1.1, xjust = 0,
       y.intersp = 1.2,
       legend = names(color_amr), fill = color_amr, title = "AMR class", 
       ncol = 1, bty = "n", title.adj = 0, cex = 2)

legend(
  x = -.05, y = 0, xjust = 0, cex = 2,
  legend = c("mOTU","acquired","functional"),
  pch = c(1,11,11),
  col = c("#838B83", "#8B0000", "#008B45"),
  pt.cex = 2, bty = "n",
  title = "Source", horiz = T,
  lwd = 2, title.adj = 0
)

legend(
  x = 0.465, y = 0, xjust = 0, cex = 2,
  legend = c("Positive", "Negative"),
  lty = c(1, 1),       
  col = c("blue", "red"),  
  lwd = 5,             
  bty = "n",           
  horiz = TRUE,
  title = "Link Weights", title.adj = 0, 
)

legend(
  x = 0.85, y = 0, xjust = 0, cex = 2,
  x.intersp = c(1,1.75),
  y.intersp = 1.75,
  legend = c("", "", ""),
  pt.cex = c(4,6,8),        
  pch = 21,                     
  pt.bg = "gray15",     
  bty = "n",           
  horiz = TRUE,
  title = "Mean CLR", title.adj = 0, 
)
dev.off()

null device 
          1 

df_bar <- gather_taxa(mg) %>%
  mutate(who = case_when(
    type == "acquired" ~ "Acquired ARGs",
    type == "functional" ~ "FG ARGs",
    #comm_id == "4" & species_uhgg == "Klebsiella grimontii" & is_from_gut == TRUE ~ "human enterobacteriaceae",
    is_from_gut ~ "Human-Gut mOTUs",
    family_ncbi == "Enterobacteriaceae" ~ "Enterobacteriaceae",
    family_ncbi == "Pseudomonadaceae" ~ "Pseudomonadaceae",
    is.na(family_ncbi) ~ "Unclassified mOTUs",
    TRUE ~ "Other mOTUs"
  )) %>%
  mutate(who = factor(who, levels = c("Acquired ARGs", "FG ARGs","Human-Gut mOTUs", 
                                      "Enterobacteriaceae", "Pseudomonadaceae",
                                      "Other mOTUs", "Unclassified mOTUs"))) %>%
  group_by(comm_id, who, color_comm) %>%
  summarise(n = n(), .groups = "drop") 

p_bar <- df_bar %>%
  mutate(comm_id = factor(comm_id, levels = as.character(1:7))) %>%
  ggplot(aes(x = comm_id, y = n, fill = who)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  scale_fill_manual(values = c("Acquired ARGs" = "#8B0000", "FG ARGs" = "#008B45",
                               "Human-Gut mOTUs" = "#CD6839", 
                               "Unclassified mOTUs" = "#D3D3D3",
                               #"human enterobacteriaceae" = "black",
                               "Pseudomonadaceae" = "#8FA4DAFF",
                               "Enterobacteriaceae" = "#E393D9FF",
                               "Other mOTUs" = "#FFFACD")) +
  theme(axis.title.y = element_blank()) +
  labs(fill = "Source") +
  xlab("Communities") +
  theme(legend.position = "right") +
  theme(
    axis.text = element_text(size = 20),
    axis.text.x = element_text(size = 20, face = "bold"),
    axis.title.x = element_text(size = 20),
    legend.text = element_text(size = 20),       
    legend.title = element_text(size = 24),     
    legend.key.size = unit(2.5, "lines"),
    plot.margin = margin(t = 100, r = 20, b = 40, l = 60)
  ) 

Merge The Plots

library(png)
library(grid)
library(ggpubr)
p_taxa <- rasterGrob(readPNG("intermediate/network_taxonomy_only.png"), interpolate = TRUE)
p_comm <- rasterGrob(readPNG("intermediate/network_community.png"), interpolate = TRUE)
p_leg <- rasterGrob(readPNG("intermediate/netw_legend.png"), interpolate = TRUE)

all_plot <- ggpubr::ggarrange(
  p_taxa,
  p_comm,
  p_leg,
  p_bar,
  
  heights = c(2, 1),
  labels = c("a", "b", "","c"),
  font.label = list(size = 32),
  label.x = c(.15,.15,0,.05),
  label.y = c(.9,.9,0,.95)
  
)
all_plot

png(filename = "../plots/networks.png", width = 2*8400, height = 2*6300, res = 600)
all_plot
dev.off()

null device 
          1 

pdf(file = "../plots/networks.pdf", width = 28, height = 21)
all_plot
dev.off()

null device 
          1 

saveRDS(mg, "cache/gs_netw_filt_with_graphics.rds")